

program define bsurvci
		syntax [if] [in] , GENerate(name) ID(string) [AT(string) AT1(string) AT2(string) AT3(string) AT4(string) TVC(string) TEXP(string) REPLACE TIES(name) strata(varname) shared(varname)  graph plotopts(string) Reps(real 1000) level(real 95)]
			
/*			if "`at'"=="" {
				if "`at1'"=="" {
					display in red "at(), or alternatively at1(), are required"
				}
			}	
*/			
			if "`strata'"=="" {
				if "`tvc'"=="" {
					if "`at'"!="" {
						bootstrap_ci_notvc `if' `in' , generate(`generate') at(`at') id(`id') `replace' ties(`ties') shared(`shared') `graph' plotopts(`plotopts') reps(`reps') level(`level')
					}
					else {
						bootstrap_ci_notvc_multiple `if' `in' , generate(`generate') at1(`at1') at2(`at2') at3(`at3') at4(`at4') id(`id') `replace' ties(`ties') shared(`shared') `graph' plotopts(`plotopts') reps(`reps') level(`level')
					}
				}
				else {
					if "`at'"!="" {
						bootstrap_ci_tvc `if' `in' , generate(`generate') at(`at') tvc(`tvc') texp(`texp') id(`id') `replace' ties(`ties') shared(`shared') `graph' plotopts(`plotopts') reps(`reps') level(`level')
					}
					else {
						display in red "Only at() allowed with tvc()"
					}
				}
			}
			else {
				if "`tvc'"=="" {
					if "`at'"!="" {
						bootstrap_ci_notvc_strata `if' `in' , generate(`generate') at(`at') id(`id') `replace' ties(`ties') strata(`strata') shared(`shared') `graph' plotopts(`plotopts') reps(`reps') level(`level')
					}
					else {
						display in red "Only at() allowed with strata()"
					}
				}
				else {
					if "`at'"!="" {
						bootstrap_ci_tvc_strata `if' `in' , generate(`generate') at(`at') tvc(`tvc') texp(`texp') id(`id') `replace' ties(`ties') strata(`strata') shared(`shared') `graph' plotopts(`plotopts') reps(`reps') level(`level')
					}
					else {
						display in red "Only at() allowed with tvc() and strata()"
					}
				}
			}
			
end
		


program define bootstrap_ci_notvc_multiple
		syntax [if] [in] , GENerate(name) AT1(string) ID(string) [REPLACE AT2(string) AT3(string) AT4(string) TIES(name) shared(varname) strata(varname) graph plotopts(string) Reps(real 1000) level(real 95)]
		
		*identify how many at* scenarios where specified
		if "`at2'"!="" {
			local at_count=2
			if "`at3'"!="" {
				local at_count=3
				if "`at4'"!="" {
					local at_count=4
				}
			}
		}
		else {
			local at_count=1
		}
		
		*for each at* scenario
		forvalues which_at=1/`at_count' {
		
			*drop variables if replace is specified
			if "`replace'"=="replace" {
				capture drop _tscurve
				capture confirm new var `generate'`which_at'
				if _rc==0 {
					display "(note: variable `generate'`which_at' not found)"
				}
				else {
					drop `generate'`which_at'
				}
				capture confirm new variable `generate'`which_at'_lb
				if _rc==0 {
					display "(note: variable `generate'`which_at'_lb not found)"
				}
				else {
					drop `generate'`which_at'_lb
				}
				capture confirm new variable `generate'`which_at'_ub
				if _rc==0 {
					display "(note: variable `generate'`which_at'_ub not found)"
				}
				else {
					drop `generate'`which_at'_ub
				}
			}

			*verify if variables already exist
			confirm new variable `generate'`which_at'_lb `generate'`which_at'_ub
		}
		
	
		*define temporal variables and names
		tempvar s0
		tempname xb

		*generate temporal file to store bootstrapped results
		tempfile data1
		tempvar `generate'_b
		preserve 
		clear
		*for each at* scenario
		forvalues which_at=1/`at_count' {
			gen ___`generate'`which_at'_b=.
		}
		gen _tscurve=.
		quietly save `data1'
		restore
		
		*determine shared() option
		if "`shared'"!="" {
			local sharedopt "shared(`shared')"
		}
		else {
			local sharedopt ""
		}

		*extract predictor variables from at1()
		local i=1
		local est_vars " "
		tokenize `at1'
		while "``i''" != "" {
			local est_vars "`est_vars' ``i''"
			local i=`i'+2
		}
		
		*calculate model and baseline survival curve for original sample
		preserve
		display as text _newline(1) "The estimation is based on the following Cox Proportional Hazards Model:"
		stcox `est_vars' `if' `in' , nohr `ties' `sharedopt' 
		matrix b=e(b)
		quietly predict `s0', basesurv

		*for each at* scenario
		local varnames 
		forvalues which_at=1/`at_count' {

			*calculate linear prediction 
			local i=1
			local j=2
			scalar `xb'=0
			local at="`at`which_at''"
			tokenize `at'
			while "``i''" != "" {
				local colnum=colnumb(b,"``i''")
				scalar `xb' = `xb' + b[1,`colnum'] * ``j''
				local i=`i'+2
				local j=`j'+2
			}
			
			*calculate survival function for scenario		
			quietly gen `generate'`which_at'=`s0'^exp(`xb')
			
			*store variable names for collapse procedure below
			local varnames `varnames' `generate'`which_at'
		}
		
		*save data for later use
		tempfile data2
		quietly gen double _tscurve=_t  `if' `in'
		sort _tscurve
		collapse (mean) `varnames', by(_tscurve)
		forvalues which_at=1/`at_count' {
			label var `generate'`which_at' ""
		}
		quietly keep if _tscurve!=.
		quietly save `data2'
		restore
		
		*run replications
		forvalues j=1/`reps' {
			*quietly {
				preserve
				bsample `if' `in' , cluster(`id')
				capture {
					stcox `est_vars' , nohr `ties' `sharedopt' 
					matrix b=e(b)
					predict `s0', basesurv
					*for each at* scenario
					forvalues which_at=1/`at_count' {
						local k=1
						local l=2
						scalar `xb'=0
						local at="`at`which_at''"
						tokenize `at'
						while "``k''" != "" {
							local colnum=colnumb(b,"``k''")
							scalar `xb' = `xb' + b[1,`colnum'] * ``l''
							local k=`k'+2
							local l=`l'+2
						}
						gen ___`generate'`which_at'_b=`s0'^exp(`xb')
					}
					rename _t _tscurve
					keep ___`generate'*_b _tscurve
					append using `data1'
					save `data1', replace
				}
				restore
			*}
			*document completed number of replications
			if `j'==1 {
			di _newline(2) "Bootstrap replications (`reps')" _newline(1) "----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5"
			}
			if _rc!=0 {
				di in red "x" _cont in smcl
			}
			else {
				di _asis "." _cont in smcl 
			}
			if mod(`j',50)==0 {
				di _asis "`j'" _newline(1) _cont
			}
		}
		
		*percentile method
		*calculate pointwise confidence intervals for each time-point
		preserve
		use `data1', clear
		quietly keep ___`generate'*_b _tscurve
		quietly keep if _tscurve!=.
		local p_low=(100-`level')/2
		local p_high=100-`p_low'
		forvalues which_at=1/`at_count' {
			bysort _tscurve: egen `generate'`which_at'_lb=pctile(___`generate'`which_at'_b), p(`p_low')
			bysort _tscurve: egen `generate'`which_at'_ub=pctile(___`generate'`which_at'_b), p(`p_high')
		}
		collapse (min) `generate'*_lb (max) `generate'*_ub, by(_tscurve)
		sort _tscurve
		forvalues which_at=1/`at_count' {
			label var `generate'`which_at'_lb "lower bound"
			label var `generate'`which_at'_ub "upper bound"
		}
		quietly save `data1', replace
		restore
		
		*merge with data from original sample
		tempvar merge
		preserve
		use `data2', clear
		quietly merge 1:1 _tscurve using `data1', gen(`merge')
		quietly sum `merge'
		if r(mean)!=3 | r(sd)!=0 {
			di _newline(2) "Note: Not all failure or censoring times were sampled during the bootstrap replications. Consider running more replications."
		}
		tempvar drop
		forvalues which_at=1/`at_count' {
			quietly gen `drop'=1 if `generate'`which_at'[_n-1]==`generate'`which_at'[_n]
			quietly replace `generate'`which_at'=. if `drop'==1
			quietly replace `generate'`which_at'_ub=. if `drop'==1
			quietly replace `generate'`which_at'_lb=. if `drop'==1
			quietly drop `drop'
		}
		rename _tscurve _t
		sort _t
		quietly save `data1', replace
		restore
		
		*merge with existing data
		tempvar original_sort
		gen `original_sort'=_n
		sort _t
		tempvar merge2
		quietly merge m:1 _t using `data1', gen(`merge2')
		sort `original_sort'
		
		*plot
		if "`graph'"=="graph" {
			local plot_list ""
			forvalues which_at=1/`at_count' { 
				quietly twoway line `generate'`which_at'_lb `generate'`which_at'_ub `generate'`which_at' _t, sort ytitle("") title("at`which_at'") lc() c(J J J) `plotopts' saving(___plot`which_at', replace) nodraw
				local plot_list "`plot_list' ___plot`which_at'.gph"
			}
			graph combine `plot_list', ycommon
			forvalues which_at=1/`at_count' { 
				erase ___plot`which_at'.gph
			}

		}


end




program define bootstrap_ci_notvc_strata
		syntax [if] [in] , GENerate(name) AT(string) ID(string) strata(varname) [REPLACE TIES(name) shared(varname) graph plotopts(string) Reps(real 1000) level(real 95)]
		
		*identify strata levels
		quietly levelsof `strata', local(stratavalues)

		
		*drop variables if replace is specified
		if "`replace'"=="replace" {
			capture drop _tscurve
			foreach s of local stratavalues { 
				capture confirm new var `generate'`s'
				if _rc==0 {
					display "(note: variable `generate' not found)"
				}
				else {
					drop `generate'`s'
				}
				capture confirm new variable `generate'`s'_lb
				if _rc==0 {
					display "(note: variable `generate'`s'_lb not found)"
				}
				else {
					drop `generate'`s'_lb
				}
				capture confirm new variable `generate'`s'_ub
				if _rc==0 {
					display "(note: variable `generate'`s'_ub not found)"
				}
				else {
					drop `generate'`s'_ub
				}
			}
		}
		
		*verify if variables already exist
		foreach s of local stratavalues { 
			capture confirm new variable `generate'`s' `generate'`s'_lb `generate'`s'_ub
		}
		
		*define temporal variables and names
		tempvar s0
		tempname xb

		*generate temporal file to store bootstrapped results
		tempfile data1
		tempvar `generate'_b
		preserve 
		clear
		foreach s of local stratavalues { 
			gen ___`generate'_b`s'=.
		}
		gen _tscurve=.
		quietly save `data1'
		restore
		
		*determine shared() option
		if "`shared'"!="" {
			local sharedopt "shared(`shared')"
		}
		else {
			local sharedopt ""
		}

		*extract predictor variables from at()
		local i=1
		local est_vars " "
		tokenize `at'
		while "``i''" != "" {
			local est_vars "`est_vars' ``i''"
			local i=`i'+2
		}
		
		*calculate model and baseline survival curve for original sample
		preserve
		display as text _newline(1) "The estimation is based on the following Cox Proportional Hazards Model:"
		stcox `est_vars' `if' `in' , nohr `ties' `sharedopt' strata(`strata')
		matrix b=e(b)
		quietly predict `s0', basesurv

		*calculate linear prediction 
		local i=1
		local j=2
		scalar `xb'=0
		tokenize `at'
		while "``i''" != "" {
			local colnum=colnumb(b,"``i''")
			scalar `xb' = `xb' + b[1,`colnum'] * ``j''
			local i=`i'+2
			local j=`j'+2
		}
		
		*calculate survival function for scenario		
		tempvar _generate
		quietly gen `_generate'=`s0'^exp(`xb')
		
		*save data for later use
		tempfile data2
		quietly gen double _tscurve=_t  `if' `in'
		sort _tscurve
		collapse (mean) `_generate', by(_tscurve `strata')
		quietly reshape wide `_generate', i(_tscurve) j(`strata')
		foreach s of local stratavalues { 
			rename `_generate'`s' `generate'`s'
			label var `generate'`s' "`strata' = `s'"
		}
		quietly keep if _tscurve!=.
		quietly save `data2'
		restore
		
		*run replications
		forvalues j=1/`reps' {
			*quietly {
				preserve
				bsample `if' `in' , cluster(`id') strata(`strata')
				capture {
					stcox `est_vars' , nohr `ties' `sharedopt' strata(`strata')
					matrix b=e(b)
					predict `s0', basesurv
					local k=1
					local l=2
					scalar `xb'=0
					tokenize `at'
					while "``k''" != "" {
						local colnum=colnumb(b,"``k''")
						scalar `xb' = `xb' + b[1,`colnum'] * ``l''
						local k=`k'+2
						local l=`l'+2
					}
					gen ___`generate'_b=`s0'^exp(`xb')
					rename _t _tscurve
					local keep_vars_b ""
					quietly levelsof `strata', local(stratavalues)
					foreach s of local stratavalues { 
						local keep_vars_b "`keep_vars_b' ___`generate'_b`s'"
					}
					collapse (mean) ___`generate'_b, by(_tscurve `strata')
					reshape wide ___`generate'_b, i(_tscurve) j(`strata')
					keep `keep_vars_b' _tscurve
					append using `data1'
					save `data1', replace
				}
				restore
			*}
			*document completed number of replications
			if `j'==1 {
			di _newline(2) "Bootstrap replications (`reps')" _newline(1) "----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5"
			}
			if _rc!=0 {
				di in red "x" _cont in smcl
			}
			else {
				di _asis "." _cont in smcl 
			}
			if mod(`j',50)==0 {
				di _asis "`j'" _newline(1) _cont
			}
		}
		
		*percentile method
		*calculate pointwise confidence intervals for each time-point
		preserve
		quietly use `data1', clear
		local keep_vars_b ""
		*quietly levelsof `strata', local(stratavalues)
		foreach s of local stratavalues { 
			local keep_vars_b "`keep_vars_b' ___`generate'_b`s'"
		}
		keep `keep_vars_b' _tscurve
		quietly keep if _tscurve!=.
		local p_low=(100-`level')/2
		local p_high=100-`p_low'
		foreach s of local stratavalues { 
			quietly bysort _tscurve: egen `generate'`s'_lb=pctile(___`generate'_b`s'), p(`p_low')
			quietly bysort _tscurve: egen `generate'`s'_ub=pctile(___`generate'_b`s'), p(`p_high')
		}
		local collapse_lb ""
		local collapse_ub ""
		foreach s of local stratavalues { 
			local collapse_lb "`collapse_lb' `generate'`s'_lb"
			local collapse_ub "`collapse_ub' `generate'`s'_ub"
		}
		collapse (min) `collapse_lb' (max) `collapse_ub', by(_tscurve)
		sort _tscurve
		foreach s of local stratavalues { 
			label var `generate'`s'_lb "`strata' = `s' lower bound"
			label var `generate'`s'_ub "`strata' = `s' upper bound"
		}
		quietly save `data1', replace
		restore
		
		*merge with data from original sample
		tempvar merge
		preserve
		quietly use `data2', clear
		quietly merge 1:1 _tscurve using `data1', gen(`merge')
		quietly sum `merge'
		if r(mean)!=3 | r(sd)!=0 {
			di _newline(2) "Note: Not all failure or censoring times were sampled during the bootstrap replications. Consider running more replications."
		}
		tempvar drop
		quietly gen `drop'=.
		foreach s of local stratavalues { 
			drop `drop'
			quietly replace `generate'`s'=`generate'`s'[_n-1] if `generate'`s'==.
			quietly gen `drop'=1 if `generate'`s'[_n-1]==`generate'`s'[_n]
			quietly replace `generate'`s'=. if `drop'==1
			quietly replace `generate'`s'_ub=. if `drop'==1
			quietly replace `generate'`s'_lb=. if `drop'==1
		}
		foreach s of local stratavalues { 
			label var `generate'`s' "`strata' = `s'"
		}
		rename _tscurve _t
		sort _t
		quietly save `data1', replace
		restore
		
		*merge with existing data
		tempvar original_sort
		gen `original_sort'=_n
		sort _t
		tempvar merge2
		quietly merge m:1 _t using `data1', gen(`merge2')
		sort `original_sort'
		
		*plot
		if "`graph'"=="graph" {
			local plot_list ""
			foreach s of local stratavalues { 
				quietly twoway line `generate'`s'_lb `generate'`s'_ub `generate'`s' _t, sort ytitle("") title("`strata' = `s'") lc(gs10 gs10) c(J J J) `plotopts' saving(___plot`s', replace) nodraw
				local plot_list "`plot_list' ___plot`s'.gph"
			}
			graph combine `plot_list', ycommon
			foreach s of local stratavalues { 
				erase ___plot`s'.gph
			}

		}


end







program define bootstrap_ci_notvc
		syntax [if] [in] , GENerate(name) AT(string) ID(string) [REPLACE TIES(name) shared(varname) graph plotopts(string) Reps(real 1000) level(real 95)]
		
		
		*drop variables if replace is specified
		if "`replace'"=="replace" {
			capture drop _tscurve
			capture confirm new var `generate'
			if _rc==0 {
				display "(note: variable `generate' not found)"
			}
			else {
				drop `generate'
			}
			capture confirm new variable `generate'_lb
			if _rc==0 {
				display "(note: variable `generate'_lb not found)"
			}
			else {
				drop `generate'_lb
			}
			capture confirm new variable `generate'_ub
			if _rc==0 {
				display "(note: variable `generate'_ub not found)"
			}
			else {
				drop `generate'_ub
			}
		}
		
		*verify if variables already exist
		confirm new variable `generate'_lb `generate'_ub
		
		*define temporal variables and names
		tempvar s0
		tempname xb

		*generate temporal file to store bootstrapped results
		tempfile data1
		tempvar `generate'_b
		preserve 
		clear
		gen ___`generate'_b=.
		gen _tscurve=.
		quietly save `data1'
		restore
		
		*determine shared() option
		if "`shared'"!="" {
			local sharedopt "shared(`shared')"
		}
		else {
			local sharedopt ""
		}

		*extract predictor variables from at()
		local i=1
		local est_vars " "
		tokenize `at'
		while "``i''" != "" {
			local est_vars "`est_vars' ``i''"
			local i=`i'+2
		}
		
		*calculate model and baseline survival curve for original sample
		preserve
		display as text _newline(1) "The estimation is based on the following Cox Proportional Hazards Model:"
		stcox `est_vars' `if' `in' , nohr `ties' `sharedopt' 
		matrix b=e(b)
		quietly predict `s0', basesurv

		*calculate linear prediction 
		local i=1
		local j=2
		scalar `xb'=0
		tokenize `at'
		while "``i''" != "" {
			local colnum=colnumb(b,"``i''")
			scalar `xb' = `xb' + b[1,`colnum'] * ``j''
			local i=`i'+2
			local j=`j'+2
		}
		
		*calculate survival function for scenario		
		quietly gen `generate'=`s0'^exp(`xb')
		
		*save data for later use
		tempfile data2
		quietly gen double _tscurve=_t  `if' `in'
		sort _tscurve
		collapse (mean) `generate', by(_tscurve)
		label var `generate' ""
		quietly keep if _tscurve!=.
		quietly save `data2'
		restore
		
		*run replications
		forvalues j=1/`reps' {
			*quietly {
				preserve
				bsample `if' `in' , cluster(`id')
				capture {
					stcox `est_vars' , nohr `ties' `sharedopt' 
					matrix b=e(b)
					predict `s0', basesurv
					local k=1
					local l=2
					scalar `xb'=0
					tokenize `at'
					while "``k''" != "" {
						local colnum=colnumb(b,"``k''")
						scalar `xb' = `xb' + b[1,`colnum'] * ``l''
						local k=`k'+2
						local l=`l'+2
					}
					gen ___`generate'_b=`s0'^exp(`xb')
					rename _t _tscurve
					keep ___`generate'_b _tscurve
					append using `data1'
					save `data1', replace
				}
				restore
			*}
			*document completed number of replications
			if `j'==1 {
			di _newline(2) "Bootstrap replications (`reps')" _newline(1) "----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5"
			}
			if _rc!=0 {
				di in red "x" _cont in smcl
			}
			else {
				di _asis "." _cont in smcl 
			}
			if mod(`j',50)==0 {
				di _asis "`j'" _newline(1) _cont
			}
		}
		
		*percentile method
		*calculate pointwise confidence intervals for each time-point
		preserve
		quietly use `data1', clear
		quietly keep ___`generate'_b _tscurve
		quietly keep if _tscurve!=.
		local p_low=(100-`level')/2
		local p_high=100-`p_low'
		bysort _tscurve: egen `generate'_lb=pctile(___`generate'_b), p(`p_low')
		bysort _tscurve: egen `generate'_ub=pctile(___`generate'_b), p(`p_high')
		collapse (min) `generate'_lb (max) `generate'_ub, by(_tscurve)
		sort _tscurve
		label var `generate'_lb "lower bound"
		label var `generate'_ub "upper bound"
		quietly save `data1', replace
		restore
		
		*merge with data from original sample
		tempvar merge
		preserve
		quietly use `data2', clear
		quietly merge 1:1 _tscurve using `data1', gen(`merge')
		quietly sum `merge'
		if r(mean)!=3 | r(sd)!=0 {
			di _newline(2) "Note: Not all failure or censoring times were sampled during the bootstrap replications. Consider running more replications."
		}
		tempvar drop
		quietly gen `drop'=1 if `generate'[_n-1]==`generate'[_n]
		quietly replace `generate'=. if `drop'==1
		quietly replace `generate'_ub=. if `drop'==1
		quietly replace `generate'_lb=. if `drop'==1
		rename _tscurve _t
		sort _t
		quietly save `data1', replace
		restore
		
		*merge with existing data
		tempvar original_sort
		gen `original_sort'=_n
		sort _t
		tempvar merge2
		quietly merge m:1 _t using `data1', gen(`merge2')
		sort `original_sort'
		
		*plot
		if "`graph'"=="graph" {
			twoway line `generate'_lb `generate'_ub `generate' _t, sort ytitle("") title("Cox proportional hazards regression") lc(gs10 gs10) c(J J J) `plotopts'

		}


end







program define bootstrap_ci_tvc_strata
		syntax [if] [in] , GENerate(name) AT(string) TVC(string) TEXP(string) ID(string) strata(varname) [REPLACE TIES(name) shared(varname)  graph plotopts(string) Reps(real 1000) level(real 95)]

		*identify strata levels
		quietly levelsof `strata', local(stratavalues)
		
		*drop variables if replace is specified
		if "`replace'"=="replace" {
			foreach s of local stratavalues { 
				capture confirm new variable `generate'`s'_lb
				if _rc==0 {
					display "(note: variable `generate'`s'_lb not found)"
				}
				else {
					drop `generate'`s'_lb
				}
				capture confirm new variable `generate'`s'_ub
				if _rc==0 {
					display "(note: variable `generate'`s'_ub not found)"
				}
				else {
					drop `generate'`s'_ub
				}
			}
		}
		
		*verify if variables already exist
		foreach s of local stratavalues { 
			capture confirm new variable `generate'`s'_lb `generate'`s'_ub
		}
		
		*generate temporal file to store bootstrapped results
		tempfile data1
		preserve 
		clear
		foreach s of local stratavalues { 
			gen ___`generate'_b`s'=.
		}
		gen _tscurve=.
		quietly save `data1'
		restore
		
		*calculate survival curve for original sample
		scurve_tvc `if' `in' , gen(`generate') at(`at') tvc(`tvc') texp(`texp') ties(`ties') shared(`shared') strata(`strata') `replace'
		tempfile data2
		preserve
		local keep_vars ""
		foreach s of local stratavalues { 
			local keep_vars "`keep_vars' `generate'`s'"
		}
		quietly keep `keep_vars' _tscurve
		quietly keep if _tscurve!=.
		quietly save `data2'
		restore
		
		*run replications
		forvalues j=1/`reps' {
			quietly {
				preserve
				bsample `if' `in' , cluster(`id') strata(`strata')
				capture {
				scurve_tvc, gen(___`generate'_b) at(`at') tvc(`tvc') texp(`texp') ties(`ties') shared(`shared') strata(`strata')
				local keep_vars_b ""
				quietly levelsof `strata', local(stratavalues)
				foreach s of local stratavalues { 
					local keep_vars_b "`keep_vars_b' ___`generate'_b`s'"
				}
				keep `keep_vars_b' _tscurve
				append using `data1'
				save `data1', replace
				}
				restore
			}
			*document completed number of replications
			if `j'==1 {
			di _newline(2) "Bootstrap replications (`reps')" _newline(1) "----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5"
			}
			if _rc!=0 {
				di in red "x" _cont in smcl
			}
			else {
				di _asis "." _cont in smcl 
			}
			if mod(`j',50)==0 {
				di _asis "`j'" _newline(1) _cont
			}
		}
		
		*percentile method
		*calculate pointwise confidence intervals for each time-point
		quietly levelsof `strata', local(stratavalues)
		preserve
		quietly use `data1', clear
		local keep_vars_b ""
		foreach s of local stratavalues { 
			local keep_vars_b "`keep_vars_b' ___`generate'_b`s'"
		}
		quietly keep `keep_vars_b' _tscurve
		quietly keep if _tscurve!=.
		local p_low=(100-`level')/2
		local p_high=100-`p_low'
		foreach s of local stratavalues { 
			bysort _tscurve: egen `generate'`s'_lb=pctile(___`generate'_b`s'), p(`p_low')
			bysort _tscurve: egen `generate'`s'_ub=pctile(___`generate'_b`s'), p(`p_high')
		}
		local collapse_lb ""
		local collapse_ub ""
		foreach s of local stratavalues { 
			local collapse_lb "`collapse_lb' `generate'`s'_lb"
			local collapse_ub "`collapse_ub' `generate'`s'_ub"
		}
		collapse (min) `collapse_lb' (max) `collapse_ub', by(_tscurve)
		*rename _tscurve _t_bootstrap
		foreach s of local stratavalues { 
			label var `generate'`s'_lb "`strata' = `s' lower bound"
			label var `generate'`s'_ub "`strata' = `s' upper bound"
		}
		quietly save `data1', replace
		restore
		
		*merge with data from original sample
		drop _tscurve `generate'*
		tempvar merge
		preserve
		quietly use `data2', clear
		quietly merge 1:1 _tscurve using `data1', gen(`merge')
		quietly sum `merge'
		if r(mean)!=3 | r(sd)!=0 {
			di _newline(2) "Note: Not all failure or censoring times were sampled during the bootstrap replications. Consider running more replications."
		}
		tempvar drop
		quietly gen `drop'=.
		foreach s of local stratavalues { 
			drop `drop'
			quietly gen `drop'=1 if `generate'`s'[_n-1]==`generate'`s'[_n]
			quietly replace `generate'`s'=. if `drop'==1
			quietly replace `generate'`s'_ub=. if `drop'==1
			quietly replace `generate'`s'_lb=. if `drop'==1
		}
		rename _tscurve _t
		sort _t
		quietly save `data1', replace
		restore
		
		*merge with existing data
		tempvar original_sort
		gen `original_sort'=_n
		sort _t
		tempvar merge2
		quietly merge m:1 _t using `data1', gen(`merge2')
		sort `original_sort'
		
		*plot
		if "`graph'"=="graph" {
			local plot_list ""
			foreach s of local stratavalues { 
				quietly twoway line `generate'`s'_lb `generate'`s'_ub `generate'`s' _t, sort ytitle("") title("`strata' = `s'") lc() c(J J J) `plotopts' saving(___plot`s', replace) nodraw
				local plot_list "`plot_list' ___plot`s'.gph"
			}
			graph combine `plot_list', ycommon
			foreach s of local stratavalues { 
				erase ___plot`s'.gph
			}

		}


end

















program define bootstrap_ci_tvc
		syntax [if] [in] , GENerate(name) AT(string) TVC(string) TEXP(string) ID(string) [REPLACE TIES(name) shared(varname) graph plotopts(string) Reps(real 1000) level(real 95)]

		*drop variables if replace is specified
		if "`replace'"=="replace" {
			capture confirm new variable `generate'_lb
			if _rc==0 {
				display "(note: variable `generate'_lb not found)"
			}
			else {
				drop `generate'_lb
			}
			capture confirm new variable `generate'_ub
			if _rc==0 {
				display "(note: variable `generate'_ub not found)"
			}
			else {
				drop `generate'_ub
			}
		}
		
		*verify if variables already exist
		confirm new variable `generate'_lb `generate'_ub

		
		*generate temporal file to store bootstrapped results
		tempfile data1
		tempvar `generate'_b
		preserve 
		clear
		gen ___`generate'_b=.
		gen _tscurve=.
		quietly save `data1'
		restore
		
		*calculate survival curve for original sample
		scurve_tvc `if' `in' , gen(`generate') at(`at') tvc(`tvc') texp(`texp') ties(`ties') shared(`shared') `replace'
		tempfile data2
		preserve
		quietly keep `generate' _tscurve
		quietly keep if _tscurve!=.
		quietly save `data2'
		restore
		
		
		*run replications
		forvalues j=1/`reps' {
			quietly {
				preserve
				bsample `if' `in' , cluster(`id')
				capture {
				scurve_tvc, gen(___`generate'_b) at(`at') tvc(`tvc') texp(`texp') ties(`ties') shared(`shared') 
				keep ___`generate'_b _tscurve
				append using `data1'
				save `data1', replace
				}
				restore
			}
			*document completed number of replications
			if `j'==1 {
			di _newline(2) "Bootstrap replications (`reps')" _newline(1) "----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5"
			}
			if _rc!=0 {
				di in red "x" _cont in smcl
			}
			else {
				di _asis "." _cont in smcl 
			}
			if mod(`j',50)==0 {
				di _asis "`j'" _newline(1) _cont
			}
		}
		
		*percentile method
		*calculate pointwise confidence intervals for each time-point
		preserve
		quietly use `data1', clear
		quietly keep ___`generate'_b _tscurve
		quietly keep if _tscurve!=.
		local p_low=(100-`level')/2
		local p_high=100-`p_low'
		bysort _tscurve: egen `generate'_lb=pctile(___`generate'_b), p(`p_low')
		bysort _tscurve: egen `generate'_ub=pctile(___`generate'_b), p(`p_high')
		collapse (min) `generate'_lb (max) `generate'_ub, by(_tscurve)
		*rename _tscurve _t_bootstrap
		label var `generate'_lb "lower bound"
		label var `generate'_ub "upper bound"
		quietly save `data1', replace
		restore
		
		*merge with data from original sample
		drop _tscurve `generate'
		tempvar merge
		preserve
		quietly use `data2', clear
		quietly merge 1:1 _tscurve using `data1', gen(`merge')
		quietly sum `merge'
		if r(mean)!=3 | r(sd)!=0 {
			di _newline(2) "Note: Not all failure or censoring times were sampled during the bootstrap replications. Consider running more replications."
		}
		tempvar drop
		quietly gen `drop'=1 if `generate'[_n-1]==`generate'[_n]
		quietly replace `generate'=. if `drop'==1
		quietly replace `generate'_ub=. if `drop'==1
		quietly replace `generate'_lb=. if `drop'==1
		rename _tscurve _t
		sort _t
		quietly save `data1', replace
		restore
		
		*merge with existing data
		tempvar original_sort
		gen `original_sort'=_n
		sort _t
		tempvar merge2
		quietly merge m:1 _t using `data1', gen(`merge2')
		sort `original_sort'
		
		*plot
		if "`graph'"=="graph" {
			twoway line `generate'_lb `generate'_ub `generate' _t, sort ytitle("") title("Cox proportional hazards regression") lc(gs10 gs10) c(J J J) `plotopts'

		}


end













program define scurve_tvc
	version 12.1
	syntax [if] [in] , GENerate(name) AT(string) TVC(string) TEXP(string) [REPLACE TIES(name) shared(varname) strata(varname) graph plotopts(string)]
	
	capture drop _tscurve


	
	if "`strata'"=="" {
		if "`replace'"=="replace" {
		capture confirm new var `generate'
		if _rc==0 {
			display "(note: variable `generate' not found)"
		}
		else {
			drop `generate'
		}
	}
	confirm new variable `generate'
	preserve
	tempvar basehc b_x H_x _xb_t
	tempfile scurve
	tempname xb coef
	capture keep `if'
	capture keep `in'
	display as text _newline(1) "Dataset has been temporarily split at failure times"
	stsplit, at(failures)
	local i=1
	local tvc_vars ""
	tokenize `tvc'
	while "``i''" != "" { 						
		capture confirm new var _``i''_t
		if _rc!=0 {
			display in red "This command needs to create interaction variables using the variables listed in at(). These new variables are called _varname_t. Please rename any existing variables which have the same name." 
			confirm new var _``i''_t
		}
		gen _``i''_t = ``i'' * `texp'
		local tvc_vars "`tvc_vars' _``i''_t"
		local ++i
	}
	local i=1
	local est_vars " "
	tokenize `at'
	while "``i''" != "" {
		local est_vars "`est_vars' ``i''"
		local i=`i'+2
	}
	if "`shared'"!="" {
		local sharedopt "shared(`shared')"
	}
	else {
		local sharedopt ""
	}
	display _newline(1) "The estimation is based on the following Cox Proportional Hazards Model:"
	stcox `est_vars' `tvc_vars', nohr `ties' `sharedopt'
	display "Note: tvc-interactions denoted by _varname_t were interacted with `texp'."
	matrix b=e(b)
	quietly predict `basehc', basehc
	sort _t
	collapse (mean) `basehc', by(_t)
	local i=1
	local j=2
	scalar `xb'=0
	tokenize `at'
	while "``i''" != "" {
		local colnum=colnumb(b,"``i''")
		scalar `xb' = `xb' + b[1,`colnum'] * ``j''
		tempname ``i''_value
		scalar ```i''_value' = ``j''
		local i=`i'+2
		local j=`j'+2
	}
	local k=1
	quietly gen `_xb_t'=0
	tokenize `tvc'
	while "``k''" != "" {
		local tvc_name "_``k''_t"
		local colnum=colnumb(b,"`tvc_name'")
		scalar `coef' = b[1,`colnum']
		quietly replace `_xb_t' = `_xb_t' + `coef' * ```k''_value' * `texp'
		local ++k
	}
	quietly gen `b_x' = 1-(1-`basehc')^exp(`xb'+`_xb_t') 
	quietly gen `H_x' = sum(`b_x')
	quietly gen `generate' = exp(-`H_x')
	keep `generate' _t
	rename _t _tscurve
	label var _tscurve "analysis time"
	label var `generate' ""
	quietly save `scurve'
	restore
	tempvar merge
	quietly merge using `scurve', _merge(`merge')
		if "`graph'"=="graph" {
		twoway line `generate' _tscurve, ytitle("") title("Cox proportional hazards regression") c(J) `plotopts'
	}
	}
	
	
	else {
	quietly levelsof `strata', local(stratavalues)
	if "`replace'"=="replace" {
		foreach s of local stratavalues { 
			capture confirm new var `generate'`s'
			if _rc==0 {
				display "(note: variable `generate'`s' not found)"
			}
			else {
				drop `generate'`s'
			}
		}
	}
	foreach s of local stratavalues { 
	confirm new variable `generate'`s'
	}
	preserve
	tempvar basehc _xb_t
	tempfile scurve
	tempname xb coef
	capture keep `if'
	capture keep `in'
	display as text _newline(1) "Dataset has been temporarily split at failure times"
	stsplit, at(failures)
	local i=1
	local tvc_vars ""
	tokenize `tvc'
	while "``i''" != "" { 						
		capture confirm new var _``i''_t
		if _rc!=0 {
			display in red "This command needs to create interaction variables using the variables listed in at(). These new variables are called _varname_t. Please rename any existing variables which have the same name." 
			confirm new var _``i''_t
		}
		gen _``i''_t = ``i'' * `texp'
		local tvc_vars "`tvc_vars' _``i''_t"
		local ++i
	}
	local i=1
	local est_vars " "
	tokenize `at'
	while "``i''" != "" {
		local est_vars "`est_vars' ``i''"
		local i=`i'+2
	}
	if "`shared'"!="" {
		local sharedopt "shared(`shared')"
	}
	else {
		local sharedopt ""
	}
	display _newline(1) "The estimation is based on the following Cox Proportional Hazards Model:"
	stcox `est_vars' `tvc_vars', nohr `ties' `sharedopt' strata(`strata')
	display "Note: tvc-interactions denoted by _varname_t were interacted with `texp'."
	matrix b=e(b)
	quietly predict `basehc', basehc
	quietly levelsof `strata', local(stratavalues)
	sort _t
	collapse (mean) `basehc', by(_t `strata')
	quietly reshape wide `basehc', i(_t) j(`strata')
	local i=1
	local j=2
	scalar `xb'=0
	tokenize `at'
	while "``i''" != "" {
		local colnum=colnumb(b,"``i''")
		scalar `xb' = `xb' + b[1,`colnum'] * ``j''
		tempname ``i''_value
		scalar ```i''_value' = ``j''
		local i=`i'+2
		local j=`j'+2
	}
	local k=1
	quietly gen `_xb_t'=0
	tokenize `tvc'
	while "``k''" != "" {
		local tvc_name "_``k''_t"
		local colnum=colnumb(b,"`tvc_name'")
		scalar `coef' = b[1,`colnum']
		quietly replace `_xb_t' = `_xb_t' + `coef' * ```k''_value' * `texp'
		local ++k
	}
	local plot_j ""
	local plot_vars ""
	foreach s of local stratavalues { 
	tempvar b_x`s' H_x`s' 
	quietly gen `b_x`s'' = 1-(1-`basehc'`s')^exp(`xb'+`_xb_t') 
	quietly gen `H_x`s'' = sum(`b_x`s'')
	quietly gen `generate'`s' = exp(-`H_x`s'')
	label var `generate'`s' "`strata' = `s'"
	local plot_j = "`plot_j' " + "J"
	local plot_vars = "`plot_vars' " + "`generate'`s'"
	}
	keep `generate'* _t
	rename _t _tscurve
	label var _tscurve "analysis time"
	quietly save `scurve'
	restore
	tempvar merge
	quietly merge using `scurve', _merge(`merge')
	if "`graph'"=="graph" {
		twoway line `plot_vars' _tscurve, title("Cox proportional hazards regression") c(`plot_j') `plotopts'
	}
	}
end



